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ABSTRACT 

We consider an important class of signal processing problems 
where the signal of interest is known to be sparse, and can be 
recovered from data given auxiliary information about how 
this data was generated. For example, a sparse green's func- 
tion may be recovered from seismic experimental data using 
sparsity optimization when the source signature is known. 
Unfortunately, in practice this information is often missing, 
and must be recovered from data along with the signal using 
deconvolution techniques. 

In this paper, we present a novel methodology to simulta- 
neously solve for the sparse signal and auxiliary parameters 
using a recently proposed variable projection technique. Our 
main contribution is to combine variable projection with spar- 
sity promoting optimization, obtaining an efficient algorithm 
for large-scale sparse deconvolution problems. We demon- 
strate the algorithm on a seismic imaging example. 

Index Terms — Sparsity optimization, variable projec- 
tion, seismic imaging 

1. INTRODUCTION 

Sparse regularization has proven to be an indispensable tool in 
many areas, including inverse problems [1 1 and compressive 
sensing |2, 3|. If a signal y is known to have a sparse or 
compressible (quickly decaying) representation y = Sx, this 
information can be used to formulate optimization problems 
of the form 

min||x|li s.t. \\ASx-b\\l<a'^, (1) 

where ^ is a measurement matrix used to measure the true 
signal y, 6 is a vector of data, and cr is a threshold that depends 
on the characteristics of measurement error In compressive 
sensing, it is possible to obtain recovery guarantees given 
properties of the true signal and A. In more general inverse 
problems, these guarantees have not been found; it is there- 
fore appropriate to consider ([T]i as a regularization approach 
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to the least squares problem. For example, in the seismic set- 
ting, where ([T]) has been particularly useful ID, A is a lin- 
earized Bom-scattering operator, S is the curvelet transform, 
and b is seismic data. While there are several popular algo- 
rithms that solve ([T), e.g. SPArsa |5|, the SPG^i ||6l algorithm 
has been particularly useful for seismic imaging fTllSllH. 

Many inverse problems contain unknown nuisance pa- 
rameters that must be estimated in order to recover the solu- 
tion [10|. In seismic imaging, the source wavelet is typically 
unknown. The main contribution of this paper is to extend the 
approach of |10 | to the sparse inversion context, and derive 
simple modifications of standard sparse solvers to incorporate 
solutions of unknown nuisance parameters on the fly. 

The paper proceeds as follows. In section |2j we intro- 
duce the seismic imaging problem with unknown wavelet, 
and formulate it as an extended sparsity promoting optimiza- 
tion problem (|3]l. In section |3] we review the ideas recently 
proposed in |10| that allow nuisance parameters to be esti- 
mated on the fly, and show how to incorporate these ideas 
into existing sparsity promoting formulations. In section |4] 
we develop an extended SPG£i algorithm to solve (|3]l, and 
we present numerical results in sectionjS] 

2. IMAGING WITH UNKNOWN WAVELET 

Seismic imaging is an approach to obtain a gridded subsur- 
face velocity perturbation y from seismic data, given a smooth 
starting model. Experiments are conducted by placing explo- 
sive sources on the surface and recording the reflected waves 
with an array of receivers on the surface. The data, di, in this 
case represents the Fourier transform of the recorded time se- 
ries for frequency i. The corresponding modeling operator, 
Fi, defines a linear relation between the recorded data for the 
i"^ frequency and the velocity perturbation. The statistical 
model for data given y is 

di=aiFiy + e„ (2) 

where is a statistical model for the measurement error, 
which is typically modeled as Gaussian, and are unknown 
complex source wavelet coefficients. Note that the model (|2]l 
is no longer linear in the decision variables {x, a) — it is bi- 
linear Since the perturbation y is known to be sparse in the 



Curvelet frame C, formulation ([TJ has been successfully used 
to recover y — Cx Il4\ when the source wavelet is known. In 
full generality, the joint inverse problem for the perturbation 
y and wavelet a is given by 



mm||a;||i s.t > 



\\d,-a,F,Cx\\l < a' 



(3) 



Note that the a parameters make the problem more difficult, 
because the forward model Q is no longer linear in the deci- 
sion variables {x, a), and the problem ^ is nonconvex. 

3. VARIABLE PROJECTION 

We begin by considering the problem 

miny^ \\di ~ aiFiCxWl s.t. ||a;||i < r . (4) 

x,a — ^ 

i 

The relationship between Q and (|3]l will be fully explained 
in section |4] In this section, we show how to use results 
from [10] to design an effective algorithm for Q. 

If we define X — {x ; < r}, problem Q is of the 

form 

V min g{x, a) , (5) 
where for any given x £ X, one can easily find 

a{x) — argmin5(a;, a) . (6) 

a 

In fact, a{x) is available in closed form when the least squares 
penalty is used in The key idea in [TOl is to consider the 
modified objective 



g{x) = g{x,a{x)), 
using the convenient formula 



(7) 



(8) 



This is basically a generalization of the variable projection 
algorithm |[lTl|. 

In the current setting this means that instead of solv- 
ing we can simply solve the modified problem 

To.i-a'^Wdi- ai{x)F^Cx\\l s.t. ||x||i < r (9) 



using e.g. the projected gradient iteration 

x''+^ ^ Px[xk - ,~g{x^)] 

with V xg computed via ([8]), with a{x) given by (|6|. By ifTOl 
Corollary 2.3], a stationary point of Q is also stationary point 
of 



4. PROJECTED REGULARIZED INVERSION 

In the previous section, we showed how to solve the extended 
problem (|4]). However, the formulation ([5]) is more impor- 
tant to us from the modeling perspective, since it is always 
easier to provide a noise threshold a than to figure out the 
'right' sparsity level r. In fact, the SPG^i algorithm solves 
formulation ([T]l by solving a series of subproblems that find 
the sparsity level r automatically given an input threshold a. 

In this section, we extend this approach to the pair of prob- 
lems (|3]l and Q. First, define 

v{T)=uiixiY\\d,-a,F,Cx\\l s.t. |la;|li<T (10) 

x^a ^ — ^ 



Suppose we find f such that w(f) = cr^. Can we expect 
that the corresponding minimizers of Q coincide with the 
minimizers of ([3])? This question is answered in surprising 
generality by lil2l Theorem 2.1]: as long as any minimizer 
a; of ([3]) satisfies \\di — aiFiCx\W = cr^, then the set of 
minimizers of (|3]l and Q match, and ||a;||i = f where v{f) = 

This general result points to using the following strategy: 
solve v{t) = by Newton's method 



v{Tk) 



(11) 



This is in fact the strategy used by SPG^i to solve the prob- 
lem ([T]l, for an appropriately defined value function. In order 
to implement this strategy, we have to be able to evaluate both 
v{t) and v'{t) for ([T0|. 

Evaluating v{t) is straightforward: we simply use the 
projected gradient method detailed in section |3] However, 
v' (r) is more difficult, since the most general variational re- 
sults for value functions lfT2l require linearity of the forward 
model, which is violated by (|2]l. Nonetheless, if we treat 
oii :— ai{x) as fixed, then by |[12i Theorem 6.2], we get 



v'{t) 



a^C^Ff{d,-a,F,Cx) 



where x solves Q for r. 

We note that the expression above is an approximation 
to the derivative, and the quality of the approximation re- 
mains to be determined. If the source weight can be esti- 
mated fairly quickly (so that it is not changing significantly 
between iterations), the approximation above becomes exact. 
For the experiments in the next section, we found that the 
proposed Newton iteration gives nearly the same result as the 
one with a fixed, 'true' source-weight. We also verified that 
when we pick cr that is reachable within our computational 
budget of 150 iterations, the algorithm correctly finds the root 
v{t) = a; see figure|4](b). 



5. NUMERICAL RESULTS 

For the experiments we use a Matlab framework for seismic 
imaging and modelling |13|, and the CurveLab toolbox lfT4l . 
Both of these are freely available for non-commercial pur- 
poses. The algorithm to solve |3] is based on the SPG^i code 
ifTSl . which is also available for download. 

We generate data for the velocity perturbation defined on 
a 201 X 301 grid with 10 m spacing depicted in figure [T] for 
6 frequencies between 5 and 25 Hz, 301 equispaced receivers 
and 15 composite sources, all located at the top of the model. 
We note that his leads to a underdetermined problem with 
27090 equations and 60501 unknowns. We use SPG£i to 
solve Q either with a fixed or with a estimated using the 
procedure outlined above. Since there is no noise in this ex- 
ample we use CT = and run the algorithm for a fixed number 
of iterations (100 in this case). 

Note that there is a fundamental non-uniqueness in the 
problem; if we multiply the source-weights with a constant 
factor, we can compensate for this by dividing the recon- 
structed model by the same factor. Therefore, we normalize 
the results such that the source-weights for each reconstruc- 
tion have the same norm (i.e., o^ii^)'^ is the same for all 
reconstructions). 

The reference result using the true source signature is 
shown in figure |2] (a). If we do not estimate the source sig- 
nature and use a; = 1, we do not get a good reconstruction, 
as is shown in figure [2] (b). Finally, if we estimate the source 
signature according to the strategy outlined in this paper, we 
obtain the result depicted in figure|2](c). 

The (normalized) optimal source weights a{x) evaluated 
at the final results as well as the true source weight are de- 
picted in figure [3] This shows that our approach is able to 
recover both the model and the source weight. The conver- 
gence histories of the SPG^i algorithm are shown in figure 

El 

6. DISCUSSION AND CONCLUSIONS 

We have proposed a novel method for estimating nuisance pa- 
rameters in the context of sparsity regularized inverse prob- 
lems, and in particular we have focused on source wavelet 
estimation in seismic imaging. The method draws on the idea 
of variable projection in order to estimate nuisance parame- 
ters on the fly, and can be implemented via a straightforward 
modification to existing sparse solvers. 

Numerical experiments demonstrate that the source wavelet 
can be recovered successfully in this manner (figure [3]), and 
that the recovery of primary parameters (specifically of the 
image) is improved when the wavelet is estimated (figure |2]). 

Note that after section [3] we can already solve Q, but 
nonetheless lot of effort is devoted in section |4] to develop a 
method for solving ([3]). The main point here is that while it is 
difficult to come up with a reasonable value for r in Q, it is 
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Fig. 1. True perturbation used for numerical experiment 

straightforward to come up with a good value for in ([3]l. In 
fact, given a finite computational budget, and no estimate for 
cr^, one can always pick ~ and perform a fixed number 
of iterations. In this mode, the algorithm in section]?] solves 
several Q problems inexactly, picking the corresponding se- 
quence of T values according to iteration ( [TT| . This is exactly 
what was done to obtain the numerical examples in section |5] 




2000, 



500 1000 1500 2000 2500 3000 
Horizontal distance (m) 



500 



£ 1000 



1500 



2000 




2000 




10 15 
Frequency (Hz) 



Fig. 3. Amplitude (a) and phase (b) of the optimal source 
weights evaluated at the models depicted in|2](b) (red) and|2] 
(c) (blue). The true source weight is also shown (dashed line). 
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Fig. 2. Reconstructed models for the true wavelet (a), a wrong 
wavelet (b) and using the wavelet estimation procedure (c). 
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Fig. 4. Convergence histories using the true wavelet (dashed), 
a wrong source weight (red) and the estimated source weight 
(blue) when (a) a — and (b) cr = 15. 
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